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Abstract 

A Monte Carlo method for finite-temperature studies of the two-dimensional 



presented. The scheme is based on an approximation which allows for an ana- 



^ ■ quantum Heisenberg antiferromagnet with random ferromagnetic bonds is 
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f~^ ■ lytic summation over the realizations of the randomness, thereby significantly 
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0^ . alleviating the "sign problem" for this frustrated spin system. The approxi- 



mation is shown to be very accurate for ferromagnetic bond concentrations of 



i-rt ■ up to ten percent. The effects of a low concentration of ferromagnetic bonds 
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O ■ on the antiferromagnetism are discussed. 
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Monte Carlo studies of frustrated quantum spin systems are difficult since positive def- 
inite weight functions cannot be constructed in general [the so called "sign problem" ] .Bel 
For random models, on the other hand, averaging over a large number of realizations of the 
randomness is necessary, which considerably increases the computational effort over what is 
required for non-random systems. Both the above difficulties are present for the antiferro- 
magnetic Heisenberg model with random ferromagnetic bonds. The two-dimensional (2D) 
version of this model is of current interest as a possible model of the copper-oxygen sheets 
of hghtly doped, but still insulating high-Tc superconductor materials. The idea, stressed by 
Aharony et al.p is that in the doped insulating phase, the holes introduced into the Cu-0 
sheets are localized at individual oxygen sites. The coupling between the copper and oxygen 
spins results in an effective ferromagnetic coupling between the copper spins adjacent to an 
oxygen spin. Due to the computational problems mentioned above, this picture has not yet 
been tested by direct numerical calculations of experimentally measurable quantities of the 
proposed model hamiltonian. 

Previous numerical work on random quantum spin systems has been largely limited to 
ID systemsQ and non-frustrated 2D models.Q Quantum Monte Carlo simulations of ran- 
dom systems with frustration have been carried out in cases where the sign problem is not 
present, such as the Ising spin glass in a transverse field.H For models with random long-range 
interactions, recent progress has been made using field-theoretic methods.l3 

In this paper a Monte Carlo method for finite-temperature studies of the 2D Heisenberg 
model with mixed antiferromagnetic and ferromagnetic nearest-neighbor couplings of equal 
strengths is presented. The scheme employs an approximation which corresponds to an 
annealing of the quenched disorder. This approximation is argued to be very accurate in the 
regime of interest for the high-T^ cuprates; a concentration of ferromagnetic bonds of less 
than ten percent. The summation over all realizations of the annealed randomness can be 
carried out analytically for each Monte Carlo configuration, thereby significantly alleviating 
the sign problem. Furthermore, in a single simulation, calculations can be carried out for 
several concentrations of ferromagnetic bonds with essentially no additional computational 



cost. 

Below, the method is described and tested for small systems. Results are presented for 
the effect of an increasing concentration of ferromagnetic bonds on the staggered structure 
factor and the uniform susceptibility. 

The model is defined by the hamiltonian 

H =Y^ JijSi ■ Sj, (1) 

where {i,j) is a pair of nearest-neighbor sites on a square lattice, and Si is a spin-| operator 
at site i. The coupling constants Jij are all of equal strength J, but their signs are random, 
with a probability p for —J (ferromagnetic) and 1 — p for +J (antiferromagnetic). In its 
current formulation the method to be presented does not allow for different ferromagnetic 
and antiferromagnetic coupling strengths. The random ±J model should, however, exhibit 
the general features associated with the presence of a low concentration of frustrating bonds. 
The computational scheme will be discussed in the context of a generalization of Hand- 
scomb's quantum Monte Carlo method,El but the same idea should be applicable to "world- 
line" methodsci as well. Consider the expectation value of an operator A at inverse temper- 
ature P = l/Zc^T: 

(A) = ^Tr{ie-^^}, Z = Tr{e-^^}. (2) 

The starting point for the generalization of Handscomb's method is to Taylor expand e~^^ 
and to write the traces in (^ as sums over diagonal matrix elements in a suitably chosen 
basis {I a)}, giving for the partition function 

^ = EE^(«I^»- (3) 

For the Heisenberg model, the basis {15*^, . . . , 5*^)}, SI G {t, i}, is chosen, and the hamil- 
tonian (|1]) is written as 

a = --Y.ab[Hi,t - H2,b] + iVfe(l - 2p)J/4, (4) 

^ 6=1 



where 

H2,b = 'S'^l(fe)'S's2{6) + 'S'si(5)5'^2{6)- (^) 

Here sl{b) and s2{b) are the sites connected by bond b, Nb = 2L^ is the number of bonds 
of the lattice, and ab is —1 if 6 is a ferromagnetic bond and +1 otherwise. The partition 
function can now be written as 

Z = 1^1^2^ ri \"|ll^a,A"/' (6) 

n=0 5„ a "'■ i=l 

where Sn denotes a sequence of n index pairs, 

^"^(::).a-fc); 

with Oj G {1, 2}, 6j G {1, . . . , A^fe} referring to an operator Ha,b- The matrix element in (^ 
is equal to or 1, and the sign of a given term is determined only by the number np of 
operators Ha,b with b being one of the ferromagnetic bonds. This sign rule is valid for a 
bipartite lattice, in which case an operator string contributing to Z must flip each spin in 
I a) an even number of times, and therefore the total number of operators H2,b must be even. 
Note that the only dependence on the realization of the randomness in (^ is in the number 
np. This is crusial in what follows. 

The actual Monte Carlo scheme has been described elsewhere,Ei and will not be discussed 
here. It suffices to note that, as has been shown above, for a given realization R of the ± J 
bonds on the lattice, an operator expectation value can be written as 

i:WcAcSc{R) 

c 
where Wc is a positive definite weight for the configuration C [C here belongs to the space 
{a^Sn-iU = 0,1,2,...} of states and index sequences], Sc{R) = (— l)"'-^'^^) is a sign which 
depends on the realization R as well as C, and Ac is a function measuring the operator A 
[the construction of Ac for various types of operators is discussed in detail in Ref. H. Here 



only operators without explicit dependence on the particular realization of the randomness 
will be considered, e.g. bulk susceptibilities and magnetic structure factors. 

In a Monte Carlo simulation the configurations C are generated using Wc as a relative 
probability distribution, and the quantities Sc{R) and AcSc{R) are measured with regular 
intervals. The expectation value of A is then given byS 

^""^^ - {Sc{R)) ' ^'^ 

and the average over the realizations of the randomness is 

where Nji is the number of realizations. If {Sc{R)) <^ 1, accurate determinations of {Sc{R)) 
and {AcSc{R)) become very time consuming. Since (S'c(-R)) in most cases approaches 
zero exponentially as the temperature is lowered, the sign problem is a severe limitation 
of the quantum Monte Carlo technique for models where one cannot construct [e.g. using 
symmetries] a weight function with a sign identically equal to one. a 

Since the weight Wc does not depend on the realization of the randomness, an estimate 
of {{A)) can be obtained by carrying out the measurements of AcSc{R) and Sc{R) on a 
set of pre-generated realizations in a single Monte Carlo simulation. Hence, the randomness 
averaging does not pose a problem. However, this by itself does not alleviate the sign 
problem, as the evaluation of the individual terms of (^^) still becomes unstable when 
Sc{R) <^ 1. The approximation introduced next will be shown to significantly reduce this 
sign problem. 

Consider the expectation values {AcSc{R)) and {Sc{R)) averaged over the randomness: 



iiSc)) = ^EiSciR)) 

^^n R 

{{AcSc)) = ^J:{AcSc{R)) (11) 

The realization-dependent averages can be written in terms of their deviations from the 
respective realization-averaged quantities as 



(AcSciR)) = {{AcSc)) + Aas{R) 

{Sc{R)) = {{Sc)) + As{R). (12) 

A realization-averaged expectation value can then be written as 

((i)) = ^^^|f|f + (^(A^), (13) 

where 0{A'^) denotes terms of order AAsiR)AsiR) and {AsiR)f. [Note that (|l|) would 
be exact if As{R) would be zero for all R, even with Aas{R) 7^ 0.] If the concentration of 



ferromagnetic bonds is low, the approximation (13) can be expected to be a good one, since 
the main contribution to ((^4)) is from realizations where the ferromagnetic bonds are far 
apart from each other. The signs {Sc{R)) should then typically be insensitive to variations 
in R. Note, however, that the approximation does contain collective impurity effects, as 
the averages in (|13|) depend in a non-trivial manner on the number of ferromagnetic bonds 
present. 

Under the above approximation, the averages over the Monte Carlo configurations and 
the realizations of the randomness have been put on an equal footing. This corresponds to 
going from quenched to annealed disorder. The computational advantage is that the averag- 
ing over the randomness can be performed analytically for each Monte Carlo configuration, 
which in effect means that each Monte Carlo measurement step corresponds to measuring 
on a very large number N^ of configurations, which for a fixed number Nf of ferromagnetic 
bonds is given by 

One might hope that this averaging enables a stable evaluation of the expectation values 
{{Sc)) and {{AcSc)) far beyond the point where estimates of {AcSc{R)) and {Sc{R)) 
become too noisy. 

Denoting by F{R) the set of ferromagnetic bonds in the realization R, the sign of a 
Monte Carlo configuration can be written as 
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sc{R)= n (-ir= n ^cw, m 

beF{R) b€F{R) 

where n^ is the number of operators acting on bond b, i.e. the number of index pairs kM 
with bi = b in Sn- Hence the sign is a product of "local signs" sc{b), with sc{b) being 
positive or negative depending on if rib is even or odd. Note that the local signs depend only 
on the Monte Carlo configuration C, and the full sign Sc{R) is calculated using only the 
local signs of the ferromagnetic bonds of R. Denoting the total number of local minus signs 
by r2„ and the total number of local plus signs by n^ = N^ — n_, the randomness averaged 
sign Sc = -^ Y.R Sc{R) of a Monte Carlo configuration is given by 



The (approximate) Monte Carlo estimate (|T3D for the disorder averaged {A) can now be 
written as 

((-4» = %^, (IT) 

where all effects of the randomness is contained in Ec, which can be easily calculated for 
each Monte Carlo configuration.E3 

Next it will be demonstrated that this estimate of ((^4)) is indeed considerably less noisy 
than (p!3D , and that the approximation involved is very good, at least when the concentration 
of ferromagnetic bonds is low. Results will be shown for the staggered structure factor 

^(^>^) = J2T.^'''^''~''''i(S^S^k)) (18) 

and the uniform susceptibility 

1 ? 

x(o, 0) = ;^ E y dT{{s;ir)sm))- (19) 

In order to test the accuracy of the "annealed" approximation, simulations oi L x L systems 
with L = 4 and 8 were carried out, and S'(7r, vr) and x(0, 0) were calculated using both (|TUp 
and ([T7|) [With (|l3), the averaging over R was done for using several hundred randomly 



generated realizations]. Fig. 1 shows results for L = 4 at temperatures T/J = 0.4,0.6, 
and 0.8. The maximum ferromagnetic bond concentration for which the averages can be 
evaluated decreases rapidly as the temperature is lowered. As expected, the approximate 
averages are easier to obtain than the exact ones. Perhaps surprisingly, no deviations of 
the approximate averages from the exact ones can be seen within statistical errors, even for 
rather high p. The antiferromagnetism is strongly suppressed by the disorder; the staggered 
structure factor decreases with p and the uniform susceptibility is enhanced. The effect 
becomes stronger as the temperature is decreased. Fig. 2 shows similar results for 8 x 
8 systems. Here the suppression of the antiferromagnetism is even stronger. Again, no 
differences between the approximate and exact results can be seen up to the maximum p for 
which they can both be relyably calculated. One might expect the errors of the annealed 
approximation to become larger at lower temperatures, where comparisons are difficult due 
to the sign problem. 

Fig. 3 shows the average sign versus the number of ferromagnetic bonds for L = 4 and 
L = 8 at T/J = 0.4, 0.6, and 0.8. There is very little size-dependence, confirming that the 
average sign at a given temperature depends essentially only on the number of ferromagnetic 
bonds present. Note that (Sc) can be accurately evaluated for the larger system even when 
it becomes extremely small.t°l Using the exact expression ([lOD is not feasible if {Sc{R)) 
becomes smaller than ^ 10"'^. This limits the maximum number of ferromagnetic bonds 
that can be studied. An accurate calculation of (Sc), on the other hand, is possible up to 
some maximum p, which is essentially independent of the system size. The evaluation of the 



expression ([T7| ) still becomes more difficult as the system size increases for operators such 
as the staggered structure factor, for which the autocorrelation time grows with the system 
size and fluctuations in (Ac^c) become problematic. 

The model (|I|) has antiferromagnetic long-range order at T = for p —>■ and fer- 
romagnetic order for p ^ 1. At intermediate p there is presumably a spin-glass phase. 
An important open question is the critical concentration of ferromagnetic bonds needed to 
destroy the antiferromagnetism. In principle flnite-size scaling of the staggered structure 
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factor can answer this question, which however is beyond the scope of this paper. Here 
some initial results for the effect of an increasing fraction of ferromagnetic bonds on systems 
of size 10 X 10 are presented. 

In Fig. 4 the staggered structure factor is graphed versus the temperature for various 
concentrations of ferromagnetic bonds. For p = 2.5% and 5%, S'(7r, tt) is significantly sup- 
pressed, but still has a temperature dependence similar to the clean system. For p = 10% 
the structure factor becomes almost temperature- independent a.t T ^ J/2. This might be 
an indication that no long-range order exists for this concentration. 

Fig. 5 shows the enhancement of the uniform susceptibility as the disorder is increased. 
For comparison, p = results for L = 64 are also shown. The finite-size effects for the 
uniform susceptibility are apparently quite small. The enhancement due to the presence of 
ferromagnetic bonds is significant already for p = 2.5%. As p is increased one would expect 
x(0, 0) to eventually diverge as T -^ 0. There are indications of such behavior for p > 5%, 
but unfortunately the sign problem limits the accuracy in this regime. 

It is unclear whether low enough temperatures can be reached for determining the critical 
concentration using the scheme presented here. It should be possible, however, to carry out 
detailed studies for larger systems at temperatures above T ^ J/4 for concentrations of a 
few percent. This should enable an assessment of the relevance of the model to the high-Tc 
cuprates. The high-temperature regime is also important in view of the recent work on 2D 
quantum antiferromagnetism based on the nonlinear sigma-model.EJ 

In summary, a scheme which alleviates the sign problem in quantum Monte Carlo sim- 
ulations of the 2D quantum Heisenberg antiferromagnet with random ferromagnetic bonds 
has been presented. Results for the staggered structure factor and the uniform susceptibility 
show that the presence of a few percent of ferromagnetic bonds substantially suppresses the 
antiferromagnetism. The method discussed here can easily be extended for 3D systems. 

I would like to thank M. Boninsegni, D. Scalapino, and M. Vekic for stimulating discus- 
sions. This work is supported by the Department of Energy under Grant No. DE-FG03- 
85ER45197. 
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FIGURES 

FIG. 1. The staggered structure factor and the uniform susceptibility versus the ferro- 
magnetic bond concentration for L = 4 at three different temperatures. Results obtained 
using ( pUp are shown as solid squares {T / J = 0.8), open circles (T/J = 0.6), and solid cir- 
cles {T/J = 0.4). The solid curves are drawn through points obtained using the "annealed" 
approximation (0). 

FIG. 2. Same as Fig. 1 for systems of size 8x8. 

FIG. 3. The average sign versus the number of ferromagnetic bonds for L = A (open 
symbols) and L = 8 (solid symbols). Circles are for T/J = 0.8, squares for T/J = 0.6, and 
triangles for T/J = OA. 

FIG. 4. The staggered structure factor for L = 10 versus the temperature for p = 
0%, 2.5%, 5%, 8%, and 10% (S decreasing with p). 

FIG. 5. The uniform susceptibility for L = 10 versus the temperature for p = 

0%, 2.5%, 5%, 8%, and 10% (x increasing with p). The dashed curve goes through p = 
results calculated for a system of size 64 x 64. 
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